Vector Autoregressive Models (VAR)

Unit Root Nonstationarity Tests

When investigating a relationship between two time series variables, we need to check whether the two models are nonstationary (ie. have unit roots). If they are, we need to decide whether or not there is a common stochastic trend. This check for stationarity will be an important step to do before any multivariate time series analysis.

Recall that a characteristic equation of an AR process must have (absolute value of) all roots > 1 in order to be stationarity. If there is a root of 1 – a “unit root” – then that process is nonstationary. This allows for unit root tests for stationarity.

Augmented Dickey-Fuller Test

The Augmented Dickey-Fuller (ADF) test is a test for a unit root in a time series. The null hypothesis is that there is a unit root (ie. nonstationarity).

The original Dickey-Fuller test was of the null hypothesis that \(\phi=1\) vs. an alternate hypothesis that \(\phi<1\) in the AR(1) model \(x_t = x_{t-1} + u_t\) (where \(u_t = w_t\) is white noise). The augmented Dickey-Fuller test expanded to allow \(u_t\) to be any stationary process, rather than strictly white noise. The ADF test approximates that stationary process with an AR model. The length of the time series determines the power of the ADF test. In R, we use the tseries::adf.test function.

According to the Cowper book, “The null hypothesis of a unit root is favoured by economists because many financial time series are better approximated by random walks than by a stationary process, at least in the short term.”

library(tseries)

## Warning: package 'tseries' was built under R version 3.4.4

data(tcm)  # built-in data on monthly treasury yields
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(tcm1y, main='Monthly yields of 1-year treasuries')  
acf(tcm1y, main=NA)
pacf(tcm1y, main=NA)

adf.test(tcm1y)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  tcm1y
## Dickey-Fuller = -2.1, Lag order = 8, p-value = 0.5
## alternative hypothesis: stationary

The plot shows this data clearly isn’t stationary. The large p-value from the ADF test tells us we do not reject the null hypothesis of unit root (ie. nonstationarity). A very persistent ACF plot, combined with a PACF plot that quickly drops off sharply are also evidence of a process having a unit root (ie. nonstationary).

Phillips-Perron Test

This procedure estimates the autocorrelations in the stationary process \(u_t\) directly, using a kernel smoother, rather than assuming the AR approximation. Because of this, it’s considered a semi-parametric test. This shares the same null hypothesis as the ADF test (that there is a unit root, ie. process is non-stationary). Critical values can be based on asymptotic theory or simulations. In R, we use tseries::pp.test.

pp.test(tcm1y)

## 
##  Phillips-Perron Unit Root Test
## 
## data:  tcm1y
## Dickey-Fuller Z(alpha) = -11, Truncation lag parameter = 6,
## p-value = 0.5
## alternative hypothesis: stationary

Again, we fail to reject the null hypothesis of a unit root (ie. non-stationarity).

Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test

In contrast to the other tests, the KPSS test has a null hypothesis of stationarity.

kpss.test(tcm1y)

## Warning in kpss.test(tcm1y): p-value smaller than printed p-value

## 
##  KPSS Test for Level Stationarity
## 
## data:  tcm1y
## KPSS Level = 3.2, Truncation lag parameter = 5, p-value = 0.01

With this test, we reject the null hypothesis of stationarity.

Multiple Time Series - Context

Classical linear regression models assume that stochastic errors are uncorrelated. That assumption is not valid for time series data. Two independent time series may actually appear to be correlated to each other, when in reality they are just coincidental, or are both driven by some third factor. This is called “spurious correlation”.

Don’t fall for temptation to fit a regression from one trending time series on another and report high \(R^2\). Common (bad) example: analysts plotting company revenue against macroeconomic time series.

The reason we do unit root tests for nonstationarity is that if two time series are independent and contain unit roots (e.g. they follow independent random walks or at otherwise nonstationary), then we may see an apparent linear relationship, due to chance similarity of the processes over the period of the time series. That apparent correlation would be spurious.

Is correlation a good measure of the dependency of two time series? No. The sample mean is only a good estimator when you can assume data is iid – not the case for time series data. Because sample variance, covariance, and correlation calculations use the sample mean, we correlation is not a good measure of dependency between the time series. High calculated “correlations” may mean nothing, or even worse, could be misleading.

Cointegration

We can also see cases wherein two individual time series do have unit roots, but are actually related. In other words, they share a common stochastic trend, rather than a spurious correlation. We refer to these as cointegrated. The formal definition of cointegration is that: two non-stationary time series \(x_t, y_t\) are cointegrated if some linear combination of them (\(a x_t + b y_t\)) is stationary.

Random Walk Example

Let’s make up a hypothetical example, with a random walk process \(\mu_t\):

\[ \mu_t = \mu_{t-1} + w_t\] Then let’s also assume we have two time series:

\[ x_t = \mu_t + w_{x,t}\]

\[ y_t = \mu_t + w_{y,t}\] where \(w_{x,t}\) and \(x_{y,t}\) are independent white noise series with zero mean. These two time series \(x_t\) and \(y_t\) are both nonstationary. However, their difference is stationary because it’s a linear combination of (stationary) white noise terms:

\[ x_t - y_t = \mu_t + w_{x,t} - \mu_t - w_{y,t} = w_{x,t} - w_{y,t}\] Hence these two time series are cointegrated because they share the same underlying stochastic trend \(\mu_t\).

Phillips-Ouliaris test for cointegration

The Phillips-Ouliaris (PO) tests the null hypothesis that the two series are not cointegrated. This is implemented in R as tseries::po.test. The function requires the series be given in matrix form.

Let’s simulate the example given above and test it:

x <- y <- mu <- rep(0, 1000) # initialize to zeros
for (t in 2:1000) {
  mu[t] = mu[t-1] + rnorm(1)
}
x <- mu + rnorm(1000)
y <- mu + rnorm(1000)

# look at a few simulated observations
head(cbind(mu,x,y))

##           mu       x       y
## [1,]  0.0000 -0.4120 -0.0281
## [2,] -0.4132 -0.9337  0.3191
## [3,]  0.4083 -1.0089  0.5152
## [4,]  0.6967  0.6669  1.3659
## [5,] -0.7741  0.1436 -0.3032
## [6,] -1.2742 -2.8530 -0.9135

# plot the time series
par(mfrow=c(1,1))
ts.plot(x)
lines(y, col='blue', lty='dashed')
lines(mu, col='red', lty='dotted')

# convert data into dataframe to do scatterplot matrix
tmp <- data.frame(x=x,y=y,mu=mu)
ggpairs(tmp)

# perform the ADF test on each:
adf.test(x)$p.value

## [1] 0.9294

adf.test(y)$p.value

## [1] 0.9433

Based on the ADF tests, we cannot reject the null hypothesis of a unit root (nonstationarity).

Let’s test stationarity of \(x_t - y_t\). It should be stationary since the difference is just a linear combination of white noise:

adf.test(x-y)

## Warning in adf.test(x - y): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  x - y
## Dickey-Fuller = -9.5, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

Based on this test, we reject the null hypothesis of a unit root, providing evidence of stationarity.

Now, check for cointegration:

po.test(cbind(x,y))  # need to bind the time series variables into a dataframe

## Warning in po.test(cbind(x, y)): p-value smaller than printed p-value

## 
##  Phillips-Ouliaris Cointegration Test
## 
## data:  cbind(x, y)
## Phillips-Ouliaris demeaned = -1000, Truncation lag parameter = 9,
## p-value = 0.01

The statistically significant p-value in the PO test tells us that we can reject the null hypothesis that they are not cointegrated, thus providing evidence for cointegration.

Vector AR Models

In VAR models, all variables are treated symmetrically, ie. they affect each other. We don’t limit ourselves to a unidirectional relationship where predictors influence the forecast variable. For example, increased income leads to increased spending, which leads to increased income, etc… In the VAR model, all variables are modeled as if they influence each other equally. If you put a lag on a particular variable in one equation, you need to put it in every equation. We call them endogenous.

All of the time series have to be stationary because it’s an extension of an AR model (VAR in level). If series are not stationary, you may need to take differences to transform it into a stationary series (VAR in difference). For a cointegrated series (where a difference of two equations is stationary), there’s something called VECM (Vector Error Correction Model). Vector ARMA model sounds theoretically possible, but in reality is near impossible to identify empirically due to the large number of parameters.

VAR Mathematical Model

Let’s look at a simple case of a VAR(1) model (\(w_{x,t}\) and \(w_{y,t}\) are bivariate white noise):

\[ x_t = \phi_{11} x_{t-1} + \phi_{12} y_{t-1} + w_{x,t}\] \[ y_t = \phi_{21} x_{t-1} + \phi_{22} y_{t-1} + w_{y,t}\] We see from the \(t\) subscripts on the white noise terms that they are contemporaneously correlated. In vector form, this can be written as:

\[ \bf{Z_t = \Phi Z_{t-1} + w_t}\] where \(\bf{Z_t} = {{x_t} \choose{y_t}}\), \(\bf{\Phi} = {{\phi_{11}\ \phi_{12}} \choose {\phi_{21}\ \phi_{22}}}\) and \(\bf{w_t} = {w_{x,t} \choose w_{y,t}}\).

This can also be expressed in terms of the backshift operator: \[[I -\Phi(B)] Z_t = w_t\]

Stationarity of a \(AR(P)\) model is defined similarly to an \(AR(p)\) model: roots of the characteristic equation must all be outside the unit circle. The characteristic equation is the determinant of the \(\bf{I -\Phi( } B \bf{)}\) matrix. This can be solved algebraically, but in R we use Mod(polyroot(...)) or abs(polyroot(...)) (equivalent).

For example, if we have the following parameter matrix:

\[ \Phi(B) =\left[ {\begin{array}{cc} 0.4 & 0.3 \\ 0.2 & 0.1 \\ \end{array} } \right] \]

The characteristic polynomial is: \[ I - \Phi(B) =\left[ {\begin{array}{cc} 1-0.4x & 0.3x \\ 0.2x & 1-0.1x \\ \end{array} } \right] \] Now, taking the determinant, we get this polynomial: \[ = (1-0.4x)(1-0.1x)-(0.3x)(0.2x) = 1 - 0.4x - 0.1x + 0.04x^2 - 0.06x^2 = 1 -0.5x - 0.02 x^2 \] Plugging in to R to solve for the roots, we get:

# these two are equivalent
Mod(polyroot(c(1,-0.5,-0.02)))

## [1]  1.861 26.861

abs(polyroot(c(1,-0.5,-0.02)))

## [1]  1.861 26.861

Since the absolute value of both roots are > 1, we know that this VAR(1) model is stationary.

Simulating and estimating a VAR(1) model

Let’s simulate two time series with these same parameters. First we have to simulate a multivariate white noise process for \(w_{x,t}\) and \(w_{y,t}\).

# first simulate multivariate white noise
library(mvtnorm)
cov_mat <- matrix(c(1, 0.8, 0.8, 1), nr=2)
w <- rmvnorm(1000, sigma=cov_mat)

# observe what the simulated multivariate white noiselooks like
head(w)

##          [,1]    [,2]
## [1,]  0.40528  0.7336
## [2,] -0.09012  0.9257
## [3,]  0.43514 -0.5123
## [4,]  0.90050  0.7476
## [5,] -1.69495 -1.9175
## [6,] -1.17387 -0.1896

tail(w)

##            [,1]    [,2]
##  [995,] -0.3269 -0.9101
##  [996,] -0.5334 -0.7130
##  [997,]  2.8682  2.8818
##  [998,] -1.2485 -1.6344
##  [999,] -0.5982 -0.8593
## [1000,]  0.2163  0.1858

plot(w)

# estimated covariance should be similar
cov(w)

##        [,1]   [,2]
## [1,] 0.8976 0.6957
## [2,] 0.6957 0.9159

# create our two multivariate white noise params
wx <- w[,1]
wy <- w[,2]
head(wx)

## [1]  0.40528 -0.09012  0.43514  0.90050 -1.69495 -1.17387

head(wy)

## [1]  0.7336  0.9257 -0.5123  0.7476 -1.9175 -0.1896

# check cross-correlation
ccf(wx, wy)

The ccf function verifies that cross-correlations are ~0 for all non-zero lags, which is required in the definition of bivariate white noise. They are contemporaneously correlated, but do not form a lag-lead relationship.

Now we’ll simulate a VAR(1) process:

x <- y <- rep(0, 1000)  # initialize to zeros
x[1] <- wx[1]
y[1] <- wy[1]
for (t in 2:1000) {
  x[t] <- 0.4 * x[t-1] + 0.3 * y[t-1] + wx[t]
  y[t] <- 0.2 * x[t-1] + 0.1 * y[t-1] + wy[t]
}

# look at ccf:
ccf(x,y, main='Cross-correlation function for a VAR(1) process')

We see the expected cross-correlation of a VAR(1) model.

# plot the scatterplot matrix
ggpairs(data.frame(cbind(x,y)))

# plot the time series
ts.plot(x)
lines(y, col='blue', lty='dashed')

# plot the ACF, PACF
par(mfrow=c(2,2))
acf(x)
acf(y)
pacf(x)
pacf(y)

We see the correlation is close to the theoretical (based on simulation parameters). Now let’s estimate a model with ar:

(xy_ar <- ar(cbind(x,y)))

## 
## Call:
## ar(x = cbind(x, y))
## 
## $ar
## , , 1
## 
##       x      y
## x 0.430 0.2355
## y 0.232 0.0472
## 
## 
## $var.pred
##       x     y
## x 0.898 0.696
## y 0.696 0.918

# xy_ar$ar[,,]

We see that the estimated parameters are reasonably close to the estimated (however, variances are very high).

Bootstrap Confidence Intervals

From Cowper:

“If the simulation is repeated many times with different realisations of the bivariate white noise, the sampling distribution of the estimators of the parameters in the model can be approximated by the histograms of the esti- mates together with the correlations between estimates. This is the principle used to construct bootstrap confidence intervals for model parameters when they have been estimated from time series.”

“The bootstrap simulation is set up using point estimates of the parameters in the model, including the variance of the white noise terms. Then time series of the same length as the historical records are simulated and the parameters estimated. A \((1−\alpha) × 100%\) confidence interval for a parameter is between the lower and upper \(\alpha/2\) quantiles of the empirical sampling distribution of its estimates.”

Example: UK and EU Exchange Rates

Load data (from Cowper) and format as time series and look at “correlation”

xrates <- read.table('us_rates.dat', header=TRUE)
head(xrates)

UK_ts <- ts(xrates$UK)
NZ_ts <- ts(xrates$NZ)
EU_ts <- ts(xrates$EU)

# plot data
ts.plot(UK_ts, NZ_ts, EU_ts, col=c('black','blue','red'))

# we're not surprised to see correlation
plot(EU_ts ~ UK_ts)

# check for unit roots
adf.test(EU_ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  EU_ts
## Dickey-Fuller = -1.6, Lag order = 10, p-value = 0.7
## alternative hypothesis: stationary

adf.test(UK_ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  UK_ts
## Dickey-Fuller = -2.3, Lag order = 10, p-value = 0.4
## alternative hypothesis: stationary

# check for cointegration
po.test(cbind(UK_ts, EU_ts))

## 
##  Phillips-Ouliaris Cointegration Test
## 
## data:  cbind(UK_ts, EU_ts)
## Phillips-Ouliaris demeaned = -22, Truncation lag parameter = 10,
## p-value = 0.04

We fail to reject the null hypothesis of the ADF tests that the series have unit roots (ie. are nonstationary). The statistically significant p-value in the PO test provides evidence of cointegration. Since they are cointegrated, we can fit a linear regression of EU on UK:

EU_UK_mod <- lm(UK_ts ~ EU_ts)
summary(EU_UK_mod)

## 
## Call:
## lm(formula = UK_ts ~ EU_ts)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.021626 -0.006835  0.000496  0.006144  0.028494 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.07437    0.00498    14.9   <2e-16 ***
## EU_ts        0.58700    0.00634    92.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.00838 on 1001 degrees of freedom
## Multiple R-squared:  0.895,  Adjusted R-squared:  0.895 
## F-statistic: 8.56e+03 on 1 and 1001 DF,  p-value: <2e-16

residualPlots(EU_UK_mod)

##            Test stat Pr(>|t|)
## EU_ts          9.822        0
## Tukey test     9.822        0

Not surprisingly, the slope coefficient is highly significant, and the \(R^2\) is very high. They are quite strongly related. We see a clear pattern in both residuals plots. Now, we’ll model the residuals as a stationary time series. Start with the typical EDA plots:

EU_UK_r <- resid(EU_UK_mod)
par(mfrow=c(2,2))
plot(EU_UK_r, main='Residuals from linear regression of UK Pound on Euro', type='l')
plot(density(EU_UK_r))  #  can plot density as an alternative to histogram
acf(EU_UK_r)
pacf(EU_UK_r)

If EU and UK data are cointegrated, we can find a linear relationship such that there’s a stationary time series. The EDA plots indicate an AR(1) model may be appropriate. Let’s use the ar function (which decides based on AIC).

(EU_UK_r_mod <- ar(EU_UK_r))

## 
## Call:
## ar(x = EU_UK_r)
## 
## Coefficients:
##      1       2       3  
##  1.043  -0.119   0.047  
## 
## Order selected 3  sigma^2 estimated as  4.08e-06

# do diagnostics on the residuals from this model
mod_resids <- na.omit(EU_UK_r_mod$resid)
plot(mod_resids)

qqnorm(mod_resids)

acf(mod_resids)

pacf(mod_resids)

The diagnostics look pretty much like white noise, although we’d rather the PACF not show any significant lags.

Example 2: Cardiovascular Mortality, Temp, and Particulate Matter

Load data from astsa package. For this example, we’ll use the vars package, too.

cmtp <- cbind(cmort, tempr, part)
head(cmtp)

##       cmort tempr  part
## [1,]  97.85 72.38 72.72
## [2,] 104.64 67.19 49.60
## [3,]  94.36 62.94 55.68
## [4,]  98.05 72.49 55.16
## [5,]  95.85 74.25 66.02
## [6,]  95.98 67.88 44.01

par(mfrow=c(1,1))
plot.ts(cmtp)  # gives a nicer display than ts.plot (which overlays them)

# look at correlation between the time series
ggpairs(data.frame(cmtp))

ccf(cmort, tempr)

ccf(cmort, part)

ccf(tempr, part)

# plot ACF and PACF plots
par(mfrow=c(3,2))
acf(cmort)
pacf(cmort)
acf(tempr)
pacf(tempr)
acf(part)
pacf(part)

par(mfrow=c(1,1))

We see they all have strong seasonality. Temperature and particulates don’t seem particularly contemporaneously correlated (based on scatterplot); the others are more so. All of the cross-correlation plots show strong lag-lead cross-correlation, too. ACF plots all show strong persistence. PACF plots also show quite a few significant lags.

Fit a VAR(1) model for the sake of a demo (not because EDA suggests it):

summary(VAR(cmtp, p=1, type='both'))  # 'both' means both constant (drift) and trend

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: cmort, tempr, part 
## Deterministic variables: both 
## Sample size: 507 
## Log Likelihood: -5116.02 
## Roots of the characteristic polynomial:
## 0.893 0.495 0.144
## Call:
## VAR(y = cmtp, p = 1, type = "both")
## 
## 
## Estimation results for equation cmort: 
## ====================================== 
## cmort = cmort.l1 + tempr.l1 + part.l1 + const + trend 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## cmort.l1  0.46482    0.03673   12.66  < 2e-16 ***
## tempr.l1 -0.36089    0.03219  -11.21  < 2e-16 ***
## part.l1   0.09942    0.01918    5.18  3.2e-07 ***
## const    73.22729    4.83400   15.15  < 2e-16 ***
## trend    -0.01446    0.00198   -7.31  1.1e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5.58 on 502 degrees of freedom
## Multiple R-Squared: 0.691,   Adjusted R-squared: 0.688 
## F-statistic:  280 on 4 and 502 DF,  p-value: <2e-16 
## 
## 
## Estimation results for equation tempr: 
## ====================================== 
## tempr = cmort.l1 + tempr.l1 + part.l1 + const + trend 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## cmort.l1 -0.24405    0.04210   -5.80  1.2e-08 ***
## tempr.l1  0.48660    0.03690   13.19  < 2e-16 ***
## part.l1  -0.12766    0.02199   -5.81  1.1e-08 ***
## const    67.58560    5.54155   12.20  < 2e-16 ***
## trend    -0.00691    0.00227   -3.05   0.0024 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 6.4 on 502 degrees of freedom
## Multiple R-Squared: 0.501,   Adjusted R-squared: 0.497 
## F-statistic:  126 on 4 and 502 DF,  p-value: <2e-16 
## 
## 
## Estimation results for equation part: 
## ===================================== 
## part = cmort.l1 + tempr.l1 + part.l1 + const + trend 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## cmort.l1 -0.12477    0.07901   -1.58     0.11    
## tempr.l1 -0.47653    0.06924   -6.88  1.8e-11 ***
## part.l1   0.58131    0.04126   14.09  < 2e-16 ***
## const    67.46350   10.39916    6.49  2.1e-10 ***
## trend    -0.00465    0.00426   -1.09     0.28    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 12 on 502 degrees of freedom
## Multiple R-Squared: 0.373,   Adjusted R-squared: 0.368 
## F-statistic: 74.7 on 4 and 502 DF,  p-value: <2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##       cmort tempr  part
## cmort 31.17  5.97  16.7
## tempr  5.97 40.96  42.3
## part  16.65 42.32 144.3
## 
## Correlation matrix of residuals:
##       cmort tempr  part
## cmort 1.000 0.167 0.248
## tempr 0.167 1.000 0.551
## part  0.248 0.551 1.000

Review the summary to see lots of interpretations about relationships between the time series, roots of characteristic equations, etc……

Now, let var package fit models with multiple values for \(p\), since we had no particular reason to pick VAR(1) before. Note that \(SC(n) = BIC\).

VARselect(cmtp, lag.max=10, type='both')

## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      9      5      2      9 
## 
## $criteria
##                1        2        3        4        5        6        7
## AIC(n)     11.74    11.30    11.27    11.23    11.18    11.15    11.15
## HQ(n)      11.79    11.38    11.38    11.37    11.35    11.35    11.38
## SC(n)      11.86    11.50    11.55    11.59    11.61    11.66    11.74
## FPE(n) 125216.92 80972.29 78268.20 75383.74 71426.10 69758.25 69749.89
##               8        9       10
## AIC(n)    11.13    11.12    11.12
## HQ(n)     11.39    11.41    11.44
## SC(n)     11.79    11.85    11.93
## FPE(n) 68122.41 67476.96 67556.45

Let’s look at the results of the \(VAR(2)\) model, recommended by BIC. Parsimonious models are usually better for forecasting.

summary(var2_fit <- VAR(cmtp, p=2, type='both'))

## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: cmort, tempr, part 
## Deterministic variables: both 
## Sample size: 506 
## Log Likelihood: -4987.186 
## Roots of the characteristic polynomial:
## 0.881 0.881 0.547 0.475 0.475 0.45
## Call:
## VAR(y = cmtp, p = 2, type = "both")
## 
## 
## Estimation results for equation cmort: 
## ====================================== 
## cmort = cmort.l1 + tempr.l1 + part.l1 + cmort.l2 + tempr.l2 + part.l2 + const + trend 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## cmort.l1  0.29706    0.04373    6.79  3.2e-11 ***
## tempr.l1 -0.19951    0.04427   -4.51  8.2e-06 ***
## part.l1   0.04252    0.02403    1.77   0.0775 .  
## cmort.l2  0.27619    0.04194    6.59  1.2e-10 ***
## tempr.l2 -0.07934    0.04468   -1.78   0.0764 .  
## part.l2   0.06808    0.02529    2.69   0.0073 ** 
## const    56.09865    5.91662    9.48  < 2e-16 ***
## trend    -0.01104    0.00199   -5.54  4.8e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 5.29 on 498 degrees of freedom
## Multiple R-Squared: 0.723,   Adjusted R-squared: 0.719 
## F-statistic:  185 on 7 and 498 DF,  p-value: <2e-16 
## 
## 
## Estimation results for equation tempr: 
## ====================================== 
## tempr = cmort.l1 + tempr.l1 + part.l1 + cmort.l2 + tempr.l2 + part.l2 + const + trend 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## cmort.l1 -0.10889    0.05067   -2.15   0.0321 *  
## tempr.l1  0.26096    0.05129    5.09  5.1e-07 ***
## part.l1  -0.05054    0.02784   -1.82   0.0701 .  
## cmort.l2 -0.04087    0.04859   -0.84   0.4007    
## tempr.l2  0.35559    0.05176    6.87  1.9e-11 ***
## part.l2  -0.09511    0.02929   -3.25   0.0012 ** 
## const    49.88049    6.85454    7.28  1.3e-12 ***
## trend    -0.00475    0.00231   -2.06   0.0399 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 6.13 on 498 degrees of freedom
## Multiple R-Squared: 0.544,   Adjusted R-squared: 0.538 
## F-statistic:   85 on 7 and 498 DF,  p-value: <2e-16 
## 
## 
## Estimation results for equation part: 
## ===================================== 
## part = cmort.l1 + tempr.l1 + part.l1 + cmort.l2 + tempr.l2 + part.l2 + const + trend 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## cmort.l1  0.07893    0.09177    0.86  0.39015    
## tempr.l1 -0.38881    0.09291   -4.18  3.4e-05 ***
## part.l1   0.38881    0.05043    7.71  6.9e-14 ***
## cmort.l2 -0.32511    0.08801   -3.69  0.00024 ***
## tempr.l2  0.05278    0.09376    0.56  0.57372    
## part.l2   0.38219    0.05306    7.20  2.2e-12 ***
## const    59.58617   12.41567    4.80  2.1e-06 ***
## trend    -0.00758    0.00418   -1.81  0.07033 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 11.1 on 498 degrees of freedom
## Multiple R-Squared: 0.468,   Adjusted R-squared: 0.46 
## F-statistic: 62.6 on 7 and 498 DF,  p-value: <2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##       cmort tempr  part
## cmort 28.03  7.08  16.3
## tempr  7.08 37.63  40.9
## part  16.33 40.88 123.4
## 
## Correlation matrix of residuals:
##       cmort tempr  part
## cmort 1.000 0.218 0.278
## tempr 0.218 1.000 0.600
## part  0.278 0.600 1.000

And do the diagnostics on the residuals:

par(mfrow=c(2,2))
ts.plot(resid(var2_fit))
qqnorm(resid(var2_fit))
acf(resid(var2_fit), lag.max=52)

pacf(resid(var2_fit), lag.max=52)

We can also test the serial correlation between the residuals with the Portmanteau Test. Its null hypothesis is that there is no autocorrelation (or serial correlation) among the residuals. The alternate hypothesis is that at least one lag has a non-zero coefficient (and thus there’s autocorrelation). The Ljung-Box test we performed on ARIMA models using box.test is also a type of Portmanteau test. However, we use serial.test for VAR models.

serial.test(var2_fit, lags.pt=12, type = 'PT.adjusted')

## 
##  Portmanteau Test (adjusted)
## 
## data:  Residuals of VAR object var2_fit
## Chi-squared = 160, df = 90, p-value = 5e-06

The highly statistically significant p-value tells us that we strongly reject the null hypothesis (of independence / no correlation between residuals).

Now, do 4-week and 8-week forecasts with this model for the sake of demonstration:

fcast <- predict(var2_fit, n.ahead=24, ci=0.95)  # 4 weeks ahead
fanchart(fcast)  # plot prediction and error

fcast <- predict(var2_fit, n.ahead=48, ci=0.95)  # 8 weeks ahead
fanchart(fcast)  # plot prediction and error